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Abstract 

We present a first-principles method for deriving effective low-energy models of electrons in 
solids having entangled band structure. The procedure starts with dividing the Hilbert space into 
two subspaces, the low-energy part ("cZ space") and the rest of the space ("r space"). The low- 
energy model is constructed for the d space by eliminating the degrees of freedom of the r space. 
The thus derived model contains the strength of electron correlation expressed by a partially 
screened Coulomb interaction, calculated in the constrained random-phase-approximation (cRPA) 
where screening channels within the d space, P^, are subtracted. One conceptual problem of this 
established downfolding method is that for entangled bands it is not clear how to cut out the d 
space and how to distinguish P^ from the total polarization. Here, we propose a simple procedure 
to overcome this difficulty In our scheme, the d subspace is cut out from the Hilbert space of the 
Kohn Sham eigenfunctions with the help of a procedure to construct a localized Wannier basis. 
The r subspace is constructed as the complementary space orthogonal to the d subspace. After this 
disentanglement, P^ becomes well defined. Using the disentangled bands, the effective parameters 
are uniquely determined in the cRPA. The method is successfully applied to 3d transition metals. 



I. INTRODUCTION 



In the last several decades many new materials with intriguing properties were discov- 
ered and synthesized. These materials range from the high-temperature superconductors to 
magnetic materials, and the latter have already found real applications in electronic indus- 
try. Typically, most of these materials contain elements from the 3d or 4/ rows and their 
electronic structure is characterized by the presence of a partially filled narrow band across 
the Fermi level. The fact that the narrow band is partially filled implies that there are many 
configurations with approximately equal weight rendering a one-particle description of the 
electronic structure problematic. Indeed, it has been recognized for a long time that many of 
the intriguing properties of these materials originate from correlations among the electrons 
residing in the partially filled narrow band. The electrons are neither fully localized, like 
core electrons, nor itinerant, like s or p electrons in alkalis or conventional semiconductors 
such as silicon or diamond. This hybrid property poses a tremendous theoretical difficulty 
for an accurate description of the electronic structure, because due to the electrons' partially 
itinerant character the problem can neither be treated as a purely atomic problem nor within 
a pure band picture. Moreover, the interaction with other electrons can be very important. 

A large amount of work has been directed to solving the correlation problem of the 
above materials. The usual approach is to consider only the narrow bands near the Fermi 
level and eliminate the degrees of freedom for the rest of the bands by the downfolding 
procedure, resulting in the well-known Hubbard model which contains an effective on-site 
Coulomb interaction, the Hubbard U. In general, the models represent multi-band systems 
containing interorbital as well as long-ranged part of Coulomb interaction. The models can 
then be solved with various sophisticated low-energy solvers such as dynamical mean-field 
theory (DMFT) [l| or solvers for lattice models 2]. 

An important issue in mapping the real system to a model Hamiltonian is how to de- 
termine the one-particle kinetic energy term and the effective interaction or the Hubbard 
U in the model. Unlike the one-particle parameters that can be downfolded from the band 
structure, the Hubbard U is much more elusive to determine and it is often treated as an ad- 
justable parameter. A widely used scheme to calculate the Hubbard U from first-principles 
is the constrained LDA (cLDA) method The cLDA method, however, is known 

from early on to yield values of U, which are too large in some cases (e.g. late transition 
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metals). It has been argued that this arises from technical difficulty in including transitions 
of electrons between the d and r space contributing the screening processes. This oversight 
leads in some cases to a larger value of U {(J. Recent extensions of the cLDA method may 
be found in Refs.j?] and 8]. Another method for determining the effective interaction is a 
scheme based on the random phase approximation (RPA). Early attempts along this direc- 



tion can be found in Refs 
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difficulty was also proposed 
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Some years ago a scheme for calculating the Hubbard U, called the constrained RPA 
(cRPA) scheme [12], was proposed. The main merit of the cRPA method over currently 
available methods is that it allows for a precise elimination of screening channels, which are 
instead to be included in a more sophisticated treatment of the model Hamiltonian. This is a 
controlled approximation without any ambiguity, expected to become asymptotically exact 
if the r-space becomes well separated from the d space. Moreover, the effective screened 
interaction can be calculated as a function of r and r', i.e., U(r, r'), independent of the basis 
functions. This allows easy access to obtaining not only on-site matrix elements but also 
off-site matrix elements as well as screened exchange matrix elements, which are usually 
taken to be the atomic value. Another merit is the possibility of obtaining the frequency- 
dependent Hubbard U, which may prove to be important. The cRPA method has now been 
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applied to a number of systems with success 0, 

Although the cRPA method is rather general, its applications to real systems have re- 
vealed a serious technical problem. The problem arises when the narrow band is entangled 
with other bands, i.e., it is not completely isolated from the rest of the bands. In many 
materials, the narrow band of interest is entangled. Even in simple materials such as the 3c? 
transition metals, the 3d bands mix with the 4s and 4p bands. Similarly, the 4/ bands of the 
4/ metals hybridize with the more extended s and p bands. For such cases, it is not clear 
anymore which part of the polarization should be eliminated when calculating the Hubbard 
U using the cRPA method. 

Some procedures to overcome the problem of determining U for entangled bands have 
been attempted. One of these is to choose a set of band indices and define the bands of 
Hubbard model as those bands corresponding to the chosen indices. Another alternative is 
to introduce an energy window and define the Hubbard bands to be those that have energy 
within the energy window. Yet another alternative is to have a combination of energy 
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window and band indices. These procedures, however, suffer from a number of difficulties. 
When choosing band indices it is inevitable that some of the states will have a character 
very different from that of the intended model. For example, in the case of 3d transition 
metals, choosing five "3d" bands will include at some k-points states which have little 3d 
character, with a considerable 4s component instead. Moreover, the chosen bands will be 
awkward to model since they do not form smoothly connected bands. A similar problem is 
encountered when choosing an energy window. A combination of band indices and energy 
window proposed in Ref. |6| partially solves the problem but it suffers from arbitrariness. 
Another procedure is, as we will discuss in detail later, to project the polarization to the 
orbitals of interest, e.g., 3d orbitals, but this procedure has been found to yield an unphysical 
result of negative static U. In this work, we offer a solution to the problem of determining 
the Hubbard U for entangled bands. The basic idea is to disentangle the narrow bands 
of interest from the rest and carry out the cRPA calculation for the disentangled band 
structure, not using the original band structure. The disentangling procedure is described 
in SecJTTl We apply the method to 3d transition metals in Sec UIII and show that the method 
is numerically stable and yields reasonable values of U. Finally the paper is summarized in 
SecflVl 

II. METHOD 

In the cRPA method we first choose a one-particle subspace {i^d}, which defines the 
model Hamiltonian, and label the rest of the Hilbert space by {i/j r }. We define to be the 
polarization within the d subspace and the total polarization is written as P. It is important 
to realize that the rest of the polarization P r = P — P^ is not the same as the polarization of 
the r subspace because it contains polarization arising from transitions between the d and r 
subspaces. Since Pd is the polarization of the model Hamiltonian, this polarization should 
be subtracted out from the total polarization when the effective parameter of the model is 
determined. The effective Coulomb interaction W r should be calculated with the rest of the 
polarization P r : 

W r (u) = [l-vP r (u)}- 1 v , (1) 
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where v is the bare Coulomb interaction. It can indeed be shown [12[| that the fully screened 
interaction is given by 



W=[l- W r P d \- x W r . (2) 

This mathematical identity ensures that W r can be interpreted as the effective interaction 
among the electrons residing in the d subspace since the screening of W r by Pd leads to the 
fully screened interaction. The matrix elements of W r in some localized functions can then 
be regarded as the frequency-dependent Hubbard U. It has been shown that the formula in 
Eq.(fTl) is formally exact, provided P r is the difference between the exact polarization P and 
the exact polarization of the d subspace P<j. In the cRPA method, P r is calculated within 
the random-phase approximation. 

If the d subspace forms an isolated set of bands, as for example in the case of the t 2g 
bands in SrV03, the cRPA method can be straightforwardly applied. However, in practical 
applications, the d subspace may not always be readily identified. An example of these is 
provided by the 3d transition metal series where the 3d bands, which are usually taken to 
be the d subspace, do not form an isolated set of bands but rather they are entangled with 
the 4s and 4p bands. To handle these cases we propose the following procedure. 

We first construct a set of localized Wannier orbitals from a given set of bands defined 
within a certain energy window. These Wannier orbitals may be generated by following the 
post-processing procedure of Souza, Marzari and Vanderbilt [if], LLtJ or other methods, such 
as the preprocessing scheme proposed by Andersen et al. within the Nth-order muffin-tin 
orbital (NMTO) method [ijj]. We then choose this set of Wannier orbitals as the generators 
of the d subspace and use them as a basis for diagonalizing the one-particle Hamiltonian, 
which is usually the Kohn-Sham Hamiltonian in the local density approximation (LDA) or 
generalized gradient approximation (GGA). The so obtained set of bands, which equivalently 
define the d subspace, may be slightly different from the original bands defined within 
the chosen energy window, because hybridization effects between the d and r spaces are 
neglected. However, it is important to confirm that the dispersions near the Fermi level well 
reproduces the original Kohn-Sham bands. From these bands we calculate the polarization 
Pd as 
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^(r)^r)^(r')^(r') _ ^(r)^(r)^rQ^*(r') j 

U — Ej + Ei + it] UJ + £j — £j — it] 

where {ipi}, {£«} (z = 1, • • • Nd) are the wavefunctions and eigenvalues obtained from diago- 
nalizing the one-particle Hamiltonian in the Wannier basis. 

It would seem sensible to define the rest of the polarization as P r = P — Pd, where P is the 
full polarization calculated using the original (Kohn-Sham) wavefunctions and eigenvalues 
{ipi}, {Ei} (i = 1,- • -N), and calculate W r according to Eq.([T]). We have found, however, 
that this procedure is numerically very unstable, resulting in some cases to unphysically 
negative static U and a large oscillation as a function of frequency. This is understandable 
given that P and Pd are obtained from two different band structures, so that low energy 
screening channels associated with the d-d transitions are not excluded from P r completely. 
Due to the singular nature of the expression in Eq. ([1]) these low-energy excitations can cause 
a large fluctuation in W r . 

Another way of calculating P r is to project the wavefunctions to the d space, 
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P d (r,r' ] u) = J2Yl 

i j 



where the projection operator V is defined as 

N d 

i=i 

The effective d polarization may be expressed as 

^*(r)^(r)^(r')^(r') ^(r)^*(r)^(r')VS*(r') 
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and P r = P — Pd can be used to calculate W r . We found that this procedure does not work 
either and is again unstable. This problem may be related to the fact that t/^'s are not 
orthogonal with each other and transitions between the states do not correspond to single 
particle-hole excitations. 

Based on these observations we propose the following procedure. We define the r subspace 

by 



\<t> i ) = (i-vM) 
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(7) 



which is orthogonal to the d subspace constructed from the Wannier orbitals. In practice it 
is convenient to orthonormalize {(pi} and prepare N — basis functions. By diagonalizing 
the Hamiltonian in this subspace a new set of wavefunctions {<pi} and eigenvalues {e^} 
(i — 1, • • • , N — Nd) are obtained. As a consequence of orthogonalizing and {ifrj}, the 
set of r bands {e{\ are completely disentangled from those of the d space {ej}, and they are 
slightly different from the original band structure {s{\. As we will see later, however, the 
numerical tests show that the disentangled band structure is close to the original one. 

The Hubbard U is calculated according to Eq.([T]) with P r = P — P^, where P is the full 
polarization calculated for the disentangled band structure. It is important to realize that 
the screening processes between the d space and the r space are included in U, although the 
d-r coupling is cut off in the construction of the wavefunctions and eigenvalues. 



III. RESULTS AND DISCUSSIONS 



As an illustration we apply the method to 3d transition metals. The electronic structure 
calculations are done in the local density approximation [3] of density functional theory 



201 ] with the full-potential LMTO implementation 21] . The wavefunctions are expanded by 



spdf + spd MTOs and a 8 x 8 x 8 k-mesh is used for the Brillouin zone summation. Spin 
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polarization is neglected. More technical details are found elsewhere 

Figure [IJa) shows the Kohn-Sham band structure of nickel. There are five orbitals having 
strong 3d character at [-5 eV:l eV], crossed by a dispersive state which is of mainly 4s 
character. Using the maximally localized Wannier function prescription with the energy 
window of [-7 eV:3 eV], interpolated "d" bands are obtained. The subsequent disentangling 
procedure gives the associated "r" bands. Comparing FigJTJb) with (a) we can see that 
there is no anti-crossing between the d bands and the r bands in (b). Otherwise the two 
band structures are nearly identical. 

In order to see the impact of the disentanglement on the screening effects, we perform 
the full RPA calculation using the disentangled band structure. The fully screened Coulomb 
interaction is compared with that for the original band structure in Figj2] where the average 
of the five diagonal terms in the Wannier basis ipi is plotted, 
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W M = \ E / dvdr'ip*{v) n {v)W{v, r'; ^)^(r')^(r') . (8) 

The two methods yield similar results The static values agree with each other within 0.2 eV, 
and the frequency dependence is weak at low frequencies. As frequency increases there is a 
sharp increase at ~20 eV, where screening by plasmons becomes ineffective. These results 
assure that screening effects can be treated accurately with the disentangled band structure. 

The Hubbard U is calculated by the constrained RPA, namely, by replacing W in Eq. dSJ) 
with W r . The results are shown in FigfJH There is no large fluctuation against frequency, 
in contrast to the methods described in SecJTTl and U(uS) shows a stable behavior. As is 
expected, U is significantly larger than W at low frequencies. This implies proper elimination 
old — d screening processes is crucial. Comparing with the previous results using a combined 
energy and band window 13J, the agreement is reasonably good. A small difference between 
these two results at low frequency may be due to a small portion of d-d screening presumably 
contained in the previous method 13), although it should be excluded from the cRPA 
calculations. 

We carried out the calculations for a series of other 3d metals as well and found in all 
the cases that (i) the present scheme is numerically stable and does not result in unphysical 
frequency dependence of U, and (ii) the value of W(u) is close to that from the original 
band structure. The latter is confirmed in FigJH where the static values {u — > limit) are 
summarized. Concerning U, the present method gives larger values compared to the previous 
results, particularly for early transition metals. Since W is nearly equal to each other in 
the two methods, the discrepancy is ascribed to the different treatment between Pd and Pd- 
We should note that Pd in the previous method depends on the choice of the window. For 
a wider window, obviously we would obtain a better agreement. Also, some states have a 
mixed character of 3d and 4s near the anti-crossing points. This makes elimination of the 
screening process difficult in the original band structure. The present scheme, on the other 
hand, enables us to determine U without ambiguity. The d bands are disentangled from the 
r bands. Consequently, the polarization in the d space is well-defined and can be removed 
completely in P r = P — Pd- 

In the present formulation, small off-diagonal matrix elements of the Kohn-Sham Hamil- 
tonian between the d space wavefunction constructed from the Wannier orbitals and the 
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r space \<f>j) are ignored. This is the reason why the anti-crossing is avoided. If the energy of 
this hybridization point in the band dispersion is smaller than the screened Coulomb energy 
and the energy scale of the interest, strictly speaking, one has to keep all of these hybridizing 
bands in the effective model, because the hybridization effects are non-perturbative. In the 
present case of transition metals, the energy crossing point of 4s and 3d bands are relatively 
larger than the screened Coulomb energy scale and the low energy models constructed only 
from the 3d Wannier orbitals may give at least qualitatively reasonable description of the 
low energy physics. 

IV. SUMMARY 

We have proposed a method to calculate the effective interaction parameters for the 
effective low-energy models of real materials when bands are entangled. The key point 
is to first properly orthogonalize the low-energy subspace contained in the models and the 
complementary high-energy subspace to each other. This orthogonalization by the projection 
technique enables the disentanglement of the bands. Once the disentangled band structure 
is obtained, the constraint RPA method can be used to determine the partially screened 
Coulomb interaction uniquely. Numerical tests for 3d metals show that the method is stable 
and yields reasonable results. The method is applicable to any system. Applications to 
more complicated systems, such as interfaces of transition metal oxides are now under way. 
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FIG. 1: (Color) (a) Kohn-Sham band structure of nickel in the LDA. (b) Disentangled band 
structure with d-r hybridization switched off. The red lines show the d states obtained by the 
maximally localized Wannier scheme, while the blue lines are disentangled r states. Energy is 
measured from the Fermi level. 
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FIG. 2: (Color online) Fully screened Coulomb interaction of nickel as a function of frequency. 
The average of the diagonal terms in the Wannier basis is plotted. The present scheme using the 
disentangled bands is compared to the results from the original band structure. 
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FIG. 3: (Color online) Hubbard U of nickel as a function of frequency. The diagonal term of 
the partially screened interaction in the Wannier basis is calculated by the present method and 
compared with the published data of Ref. [yj . 
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FIG. 4: (Color online) Static values of (a) fully screened Coulomb interaction W and (b) Hubbard 
U for 3d metals. 
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